/*  snmath.h  by K.Tsuru  */
/*********************
SN library
Mathematical functions since version 2.21
Separate from "snfunc.h" for convenience.
**********************/
#ifndef SN_MATH_H
#define SN_MATH_H

#ifndef SN_FUNCTIONS_H
#include "snfunc.h"
#endif // SN_FUNCTIONS_H

/////////////////// common class /////////////////
// x to n-th power i.e. return x^n since version 2.30
// unsigned long version
template <class T> T PowUL(const T& x, ulong n) {
    if(!n) return T(1.0);
    T y(1), z(x);

    while(1){
        if( n & 1 ) y *= z;
        n /= 2;
        if( !n ) break;
        z *= z;
    }
    return y;
}
// long version
template <class T> T PowL(const T& x, const long n) {
    if(!n) return T(1.0);
    T y(1), z(x);
    long m = labs(n);
    while(1){
        if( m & 1 ) y *= z;
        m /= 2;
        if( !m ) break;
        z *= z;
    }
    if(n > 0) return y;
    return 1/y;
}
/********************************************
 Tailor series expansion for sine and cosine  since ver. 2.31
 Type : double, long double or SDouble can be  used.
 for SDouble "eps" should be cast "SDouble(1.0e-240)"
 but do not recomend.
 Please use Sin(x) or Cos(x).
*********************************************/
template <class Type>Type sinTaylorSeries(const Type& x, const Type& eps){
    int k = 2;
    Type sum = x, xsq = x*x, term = x;
    do {    
        term *= (-1)*xsq/((k+1)*k);
        sum += term;
        k += 2;
    } while (abs(term) > eps);
    return sum;
}
template <class Type>Type cosTaylorSeries(const Type& x, const Type& eps){
    int k = 1;
    Type sum = 1.0, xsq= x*x, term = 1.0;
    do {    
        term *= (-1)*xsq/((k+1)*k);
        sum += term;
        k += 2;
    } while (abs(term) > eps);
    return sum;
}
////////////// SLong class /////////////////
// greatest common divisor
SLong gcdL(const SLong& a, const SLong& b); //  2000
inline SLong gcd(const SLong& a, const SLong& b) { return gcdL(a, b); }
// least common multiple
SLong lcmL(const SLong& a, const SLong& b); // 2008 since ver 2.8
inline SLong lcm(const SLong& a, const SLong& b){ return lcmL(a, b); }

SLong Lpow(const SLong& x, ulong n);    // n-th power of x(x^n) 2002
//10^n
SLong Lpow10(long n);                               // "friend" is deleted since version 2.192

/**********************************************
It returns a SLong random number in radix.
When size = 0 it is taken as minArraySize.
If 'size' has a value is greater than the maximum 
size it is adjusted to the maximum one.
***********************************************/
SLong LRand(uint size, fType radix = DRADIX);           //  2003
// square root of SLong integer N, q = [sqrt(N)]
// sgn has sign of N - q*q, then sgn >= 0
// If N is a square of an integer, sgn = 0.
SLong LSqrt(const SLong& N, int* sgn = NULL);       //  2004
// return M^d mod n
SLong MpowMod(const SLong& M, ulong d, const SLong& n);// 2005
inline SLong Pow2(ulong p){ return Lpow(2, p); }  // inline function since ver. 2.3
// If m and n have a common divisor in type 2^b*R^r, they are divided by it.
void  ReduceByPow2Rdx(SLong& m, SLong& n, SLong* divisor = NULL);// 2007

inline SLong Fact(ulong n) { return FactPF(n); }   // factorial n!  since version 2.21
const ulong nMax_combL = 2000UL; // maximum value of n
SLong combL(ulong n, ulong k); // combination number(binomial) nCk using SInteger arithmertic that is fast for small n, k 4100
ulong combUL(ulong n, ulong k); // ulong version(an assistant function) n < 30 (32 bits)
SLong comb(ulong n, ulong k);  // combination number(binomial) nCk  4102
inline SLong binomial(ulong n, ulong k) { return comb(n,k); }// same as above
inline SLong perm(const ulong n, const ulong r) { return permL(n, r); } // permutation number nPr "PrimeFactor" method
SLong  EulerNum(uint n);  // n-th Euler number              4101
SLong  TanCoeff(uint n);  // n-th tangent coefficient       4104
void TanCoeffFree();

void   TanCoeffTable(SNBlock <SLong>& T, uint N);   //table of tangent coefficient 4105

/////////////////// SFraction class /////////////////
SFraction BernNum(uint n);      // Bernouilli's number 7100
void      BernNumTable(SNBlock <SFraction>& Bn, uint N); // table of Bernouilli's number 7101

// Digit figures of factorial n! by Stirling's formula
unsigned long Stirling(unsigned long n);    // ver. 2.17.1 change into "unsigned long" from "long".

/////////////////// SDouble class /////////////
// Convert SDouble to an approximate in double.
// If set err = 0, return DBL_MAX for overflow or zero for underflow error in double.
// Error ID(OVERFLOW_ERR or UNDERFLOW_ERR) can be obtained by SNManeger::SNError().
double  doubleD(const SDouble& sd, int err = 1);// convert to double 3001
ldouble ldpow10(const int n); // since version 2.31
ldouble doubleLD(const SDouble& sd, int err = 1); // convert to long double since 2.31
SDouble Dpow(const SDouble& x, long n); //n-th power of x(x^n)  3002
inline SDouble Dpown(const SDouble& x, long n) { return Dpow(x, n); }  // since ver. 2.30
SDouble DpowD(const SDouble& x, const SDouble& p);  // x^p = exp(p*log(x)) for x > 0 3003
inline SDouble Dpow10(long n) {
    SDouble u(1.0); return u.MultPow10(n);
}
// inline SDouble Dpow10(const SDouble& p); // 10^p since ver 2.30 see below for the definition

/************************************************
It returns a random number with radix=rdx and exponent=exp.
Default value exp = 0, radix = DRADIX.
If size = 0 or size is greater than maxmum size it is
adujusted to the maximum size.
*************************************************/
SDouble DRand(uint size = 0, int exp = 0, fType rdx = DRADIX);// 3005  The order of argument is changed since ver. 2.30
inline SDouble Drand(uint size = 0, int exp = 0, fType rdx = DRADIX) { // small capital version
    return DRand(size, exp, rdx);
}
/*-----------------
divide x into integer and decimal part
SDouble version of modf()
x = ipart + dpart
For example
x = -1.5
ipart = -1.0
It returns  -1.5 - (-1) = -0.5.
--------------------*/
SDouble Modf(const SDouble& x, SDouble& ipart); //  3008

// SDouble version of fmod()
SDouble Fmod(const SDouble& x, const SDouble& y, SDouble& ipart); //  3009
// convert SDouble x to E form i.e. mant*10^exp
// It returns mantissa.
SDouble SDtoE_Form(const SDouble& x, long* exp);
// square root
SDouble RecSqrt(const SDouble& x);      // 1/sqrt(x)    3006
SDouble Sqrt(const SDouble& x);         // sqrt(x)      3007
// cubic root
SDouble Cbrt(const SDouble& x);         // x^(1/3)      3010
// reciporocal cubic root
SDouble RecCbrt(const SDouble& x);      // x^(-1/3)     3011
// quartic root
SDouble Qrrt(const SDouble& a);         // x^(1/4)      3012
// reciporocal quartic root
SDouble RecQrrt(const SDouble& a);      // x^(-1/4)     3013

// exponential and hyperbolic functions
SDouble Expower(long n);        // n-th power of e      3302
/**************************************************************************
trigonometric functions calculated by DRAIDX
**************************************************************************/
SDouble Acos(const SDouble& x); //                      3101
SDouble Asin(const SDouble& x); //                      3102
const uint thFigAtan = 400u; // threshold value
inline SDouble Atan(const SDouble& x) { return (x.EffFig() > thFigAtan) ? AtanA(x) : AtanT(x); } // since version 2.30
inline SDouble Cos(const SDouble& x) { return CosBS(x); } // since version 2.31
inline SDouble Sin(const SDouble& x) { return SinBS(x); }// since version 2.31
inline SDouble Tan(const SDouble& x) { return TanBS(x); }// since version 2.30
// SDouble version of C's "double atan2(double y, double x)" arctan(y/x)  3110
SDouble Atan2(const SDouble& y, const SDouble& x);
/**************************************************************************
It provides the exponential function exp(x).
From the condition "exp(x) < DRADIX^DRADIX_EXP_MAX" the range of 'x' must be
within
 |x| < DFIGURES*log(10)*DRADIX_EXP_MAX = 301759.17...
to inline since version 2.21-2.30
***************************************************************************/
const uint thFigExpx = 510u; // threshold value
inline SDouble Exp(const SDouble& x) {  // since ver 2.30
    return (x.EffFig() > thFigExpx ) ? ExpBS(x) : ExpDS(x);
}

inline SDouble Cosh(const SDouble& x) {
    return CoshBS(x);
}
inline SDouble Sinh(const SDouble& x) {
    return SinhBS(x);
}
inline SDouble Tanh(const SDouble& x){
    return TanhBS(x);
}
// inverse hyperbolic functions
SDouble Acosh(const SDouble& x);                     // 3309
SDouble Asinh(const SDouble& x);                     // 3310
SDouble Atanh(const SDouble& x);                     // 3311

// logalithm functions
inline SDouble Log(const SDouble& x) { // log(x) since ver 2.30
    return LogE(x);
}
inline SDouble Log10(const SDouble& x){ return Log(x)*RecLog10(); } // log_10(x)
// y = a^x, log y = x
inline SDouble Pow(const SDouble &a, const SDouble &x) { return DpowD(a, x); }
/***********************************************************************
Prototype declaration of constant functions has a formula
SDouble func(const SDouble* userV = NULL,SDouble (*pfCalcFunc)() = fname);
userV      : a value given by user
pfCalcFunc : default function, NULL is ok. A funcion which is faster than that
of SN library can be used.
See also "SetConstByFile()" in below.
************************************************************************/
SDouble E(const SDouble* e = NULL,
    SDouble (*pfCalcFunc)() = SNE); // base of natural logarithm 3521
void EFree();   //free the momory of exp1
uint ESize();

SDouble Log10(const SDouble* ln10 = NULL,
    SDouble (*pfCalcFunc)() = SNLog10);   // log(10.0)  3522
uint Log10Size();
void Log10Free();

SDouble Pi(const SDouble* pi = NULL,
    SDouble (*pfCalcFunc)() = SNPi);    // pi           3523
uint PiSize();
SDouble M2Pi();         // 2/pi                     3506
void M2PiFree();            //free memory of m2pi
/** 1/e If "invE=true", it returns "1/E()" else evaluates by
the series(BS method for test "RecE()*E()=1 ?"). 3507 **/
SDouble RecE(bool invE=true);
void RecEFree();        //free the memory of recE

SDouble MPi2();         // pi/2                     3511
void MPi2Free();            //free the memory of mpi2

SDouble MPi4();         // pi/4                     3550
void MPi4Free();            //free the memory of mpi4

SDouble RecLog10();     // 1/log(10.0)              3551
void recln10Free();     //free the memory of recln10

inline SDouble Dpow10(const SDouble& p) {   // 10^p since ver 2.30 It needs the definition "Log10()".
    return Exp(p * Log10());                // There is not "double pow10(double x);" in new gcc.
}

/////////////////// SComplex class /////////////////
#ifndef SCOMPLEX_H
#include "scomplex.h"
#endif // SCOMPLEX_H

inline SComplex Conj(const SComplex& z) {
    return SComplex(z.Real(), -z.Imag());
} // to inline since version 2.21
 // |z| most simple. If you need more accurate value, please improve your precision or use CabsH(z) below.
inline SDouble  Cabs(const SComplex& z) { return Sqrt(z.Norm()); }
SDouble CabsH(const SComplex& z); // |z| high precision version
inline SDouble Arg(const SComplex& z) { return Atan2(z.Imag(), z.Real()); }  // argument in xy plane
inline SDouble Real(const SComplex& z){ return z.Real();}   // real part
inline SDouble Imag(const SComplex& z){ return z.Imag();}   // imaginary part
inline SDouble Re(const SComplex& z)  { return z.Real();}   // Re(z) real part
inline SDouble Im(const SComplex& z)  { return z.Imag();}   // Im(z) imaginary part
inline SDouble Norm(const SComplex& z){ return z.Norm(); }  // square of the magnitude
inline SComplex MultI(const SComplex& z) { return SComplex(-Im(z), Re(z)); } // i*z
inline SComplex MultMI(const SComplex& z) { return SComplex( Im(z), -Re(z) ); } // -i*z = z/i
// Convert SComplex to an approximate in complex <double>.
inline complex <double> complexD(const SComplex &z) {
    return complex <double> (doubleD(z.Real()), doubleD(z.Imag()));
}
/*********************************
SComplex mathematical functions
Acooding to GCC's "complex.h" the initial letter
of funtion name is 'C';
**********************************/
SComplex Csqrt(const SComplex& z);              // sqrt(z) 9101
inline SComplex Cpolar(const SDouble& r, const SDouble& angle) {// r*{cos(angle)+i*sin(angle)}
//  return SComplex( r*Cos(angle), r*Sin(angle) );
    return CpolarBS(r, angle); 
}
SComplex Cpow(const SComplex& z, int n);        // z^n 9103
SComplex Cpow(const SComplex& x, const SDouble& y);// x^y 9104
inline SComplex Cpow(const SComplex& x, double y) {
    return Cpow(x, (SDouble)y);
}
SComplex Cpow(const SDouble& x, const SComplex& y);// x^y 9105
inline SComplex Cpow(double x, const SComplex& y) {
    return Cpow( (SDouble)x, y);
}
SComplex Cpow(const SComplex& x, const SComplex& y);// x^y 9106

SComplex Cexp(const SComplex& z);               // exp z 9107
SComplex Clog(const SComplex& z);               // log z 9108

SComplex Ccosh(const SComplex& z);              // cosh z 9109
SComplex Csinh(const SComplex& z);              // sinh z 9110
SComplex Ctanh(const SComplex& z);              // tanh z 9111

SComplex Ccos(const SComplex& z);               // cos z 9112
SComplex Csin(const SComplex& z);               // sin z 9113
SComplex Ctan(const SComplex& z);               // tan z 9114

//The following are added since ver. 2.18.
SComplex Cacos(const SComplex& z);              // arccos z 9115
SComplex Casin(const SComplex& z);              // arcsin z 9116
SComplex Catan(const SComplex& z);              // arctan z 9117

SComplex Cacosh(const SComplex& z);             // arcosh z 9118
SComplex Casinh(const SComplex& z);             // arsinh z 9119
SComplex Catanh(const SComplex& z);             // artanh z 9120

#endif // SN_MATH_H